En este documento se aplica la regresión Poisson para predecir el número de vehículos registrados diariamente en Colombia en el Registro Único de Tránsito (RUNT).
Los datos son obtenidos del RUNT a través de la compra de información y se muestran aquí solo con fines educativos.
El número de unidades registradas en el RUNT está relacionada con la actividad comercial y de las oficinas de tránsito, que operan mayoritariamente entre lunes y viernes. Además, las fechas especiales pueden tener un impacto.
La base de datos original tiene solo dos campos: fecha y unidades registradas en el RUNT en esa fecha. La base de datos se aumenta utilizando variables indicadoras para los días de la semana, los meses del año y las fechas especiales. Esta es la base que se lee a continuación.
Antes de leer los datos de un archivo externo, y si estos no ocupan un volumen importante, es conveniente visualizarlos para saber qué tipos de datos se tienen y cómo leerlos mejor.
En este caso se leerán los datos con el paquete readr y la función read_delim():
library(readr)
datos_ventas_autos <- read_delim("datos_ventas_autos.csv",
";", escape_double = FALSE, trim_ws = TRUE)
head(datos_ventas_autos) # muestra una vista de los primeros 10 registros
La función header() muestra una vista de los primeros 6 registros.
El análisis descriptivo busca visualizar comportamientos de los datos como dispersión, simetría o asimetría, relaciones entre variables, etc., mediante el cálculo de unas variables de resumen, llamadas estadísticos, y la representación gráfica de los datos.
Para guiar el análisis de la información resulta útil representar gráficamente los datos o una muestra de los mismos. En este caso, como la variable a predecir es el número de unidades registradas diariamente en el RUNT y esta variable depende del tiempo, se desea representarla como una serie temporal.
Sin embargo, la variable Fecha es una variable tipo caracter:
class(datos_ventas_autos$Fecha)
## [1] "character"
Para representar el número de unidades registradas en el RUNT en el tiempo es útil cambiar el tipo de variable a caracter. El siguiente código hace esto usando la función as.Date():
datos_ventas_autos$Fecha<-as.Date(datos_ventas_autos$Fecha,"%d/%m/%Y")
Ahora se procede a graficar los datos para darse una idea de algunos comportamientos:
library(plotly)
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~Fecha,
y = ~und,
type = "scatter" ,mode = "lines",
line=list(width=1,color='rgb(205, 12, 24)'))%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día"),
yaxis=list(title="Unidades"))
Ahora veamos el comportamiento anual. Creemos la columna “year” (año):
datos_ventas_autos$year<-format(datos_ventas_autos$Fecha,"%Y")
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~Fecha,
y = ~und,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1))%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día"),
yaxis=list(title="Unidades"))
Otra forma de analizar los resultados por año es con la función aggregate(). Utilicemos esta función para obtener el promedio diario para cada año de registros de vehículos:
aggregate(und~year,data=datos_ventas_autos,FUN=mean)
Ahora obtengamos el promedio diario para cada mes y cada año:
datos_ventas_autos$mes<-strftime(datos_ventas_autos$Fecha, format = "%B")
datos_ventas_autos$mes<-factor(datos_ventas_autos$mes,levels = month.name)
aggregate(und~year*mes,data=datos_ventas_autos,FUN=mean)
Veamos este promedio representado en un gráfico:
aggregate(und~year*mes,data=datos_ventas_autos,FUN=mean)%>%
plot_ly(x = ~mes,
y = ~und,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1))%>%
layout(title='Promedio diario mensual de registros el en RUNT',
xaxis=list(title="Día"),
yaxis=list(title="Unidades"))
Ahora obtengamos el día del año para comparar todos los años:
datos_ventas_autos$diaano<-strftime(datos_ventas_autos$Fecha, format = "%j")
Obtengamos la gráfica del registro de vehículos diarios para cada año:
library(plotly)
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~diaano,
y = ~und,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1))%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día"),
yaxis=list(title="Unidades"))
Ahora veamos el comportamiento mensual. Primero se extrae la información del mes y se expresa como una variable categórica ordenada:
datos_ventas_autos$diames<-format(datos_ventas_autos$Fecha,"%d")
La siguiente función “automatiza” los pasos para generar un gráfico para cada mes:
grafico<-function(month,data,sl=FALSE){
plot_ly (data=subset(data,subset = (Fecha<="2017-12-31" & mes==month)),
x = ~diames,
y = ~und,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1),
legendgroup=~year,
showlegend = sl)%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día"),
yaxis=list(title="Unidades"))->p
return(p)
}
Apliquemos esta función a todos los meses del año:
meses<-unique(datos_ventas_autos$mes)
gr1<-grafico(meses[1],data=datos_ventas_autos,sl=FALSE)
gr2<-grafico(meses[2],data=datos_ventas_autos,sl=FALSE)
gr3<-grafico(meses[3],data=datos_ventas_autos,sl=FALSE)
gr4<-grafico(meses[4],data=datos_ventas_autos,sl=FALSE)
gr5<-grafico(meses[5],data=datos_ventas_autos,sl=FALSE)
gr6<-grafico(meses[6],data=datos_ventas_autos,sl=FALSE)
gr7<-grafico(meses[7],data=datos_ventas_autos,sl=FALSE)
gr8<-grafico(meses[8],data=datos_ventas_autos,sl=FALSE)
gr9<-grafico(meses[9],data=datos_ventas_autos,sl=FALSE)
gr10<-grafico(meses[10],data=datos_ventas_autos,sl=FALSE)
gr11<-grafico(meses[11],data=datos_ventas_autos,sl=FALSE)
gr12<-grafico(meses[12],data=datos_ventas_autos,sl=TRUE)
subplot(gr1,gr2,gr3,gr4,gr5,gr6,gr7,gr8,gr9,gr10,gr11,gr12,nrows = 4,shareX = TRUE,
shareY = TRUE,which_layout=1)
Una alternativa, más simple pero con menos interacción, para la gráfica anterior se obtiene usando la función xyplot():
library(lattice)
xyplot(und~diames|mes,groups = year,data=datos_ventas_autos,
subset = (Fecha<="2017-12-31"),type=c("l","g"),
auto.key = list(lines = TRUE, space = "bottom",columns=7,title="Año",cex=0.5),
layout=c(3,4))
Ahora utilicemos el diagrama de caja y bigotes para explorar relaciones:
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~year,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Año"),
yaxis=list(title="Unidades"))
Ahora, hacemos algo similar a lo anterior para los meses del año y los días de la semana:
datos_ventas_autos$dia_semana<-weekdays(datos_ventas_autos$Fecha)
datos_ventas_autos$dia_semana<-ordered(datos_ventas_autos$dia_semana,levels=c( "Tuesday", "Wednesday", "Thursday",
"Friday", "Saturday", "Sunday","Monday"))
Veamos el diagrama de caja y bigotes para cada mes:
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~mes,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Mes"),
yaxis=list(title="Unidades"))
Veamos el diagrama de caja y bigotes para cada día de la semana:
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~dia_semana,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día de la semana"),
yaxis=list(title="Unidades"))
¿Cuál será el efecto de los días feriados?
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~Feriado,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Día feriado (1 es Si)"),
yaxis=list(title="Unidades"))
¿Cuál será el efecto del día del amor y la amistad?
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~AmoryAmistad,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Amor y amistad (1 es Si)"),
yaxis=list(title="Unidades"))
¿Cuál será el efecto de la Semana Santa?
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
x = ~`Semana Santa`,
y = ~und,
type = "box")%>%
layout(title='Registros el en RUNT',
xaxis=list(title="Semana santa (1 es Si)"),
yaxis=list(title="Unidades"))
Se ajustarán modelos con la información disponible hasta el 31 de diciembre de 2016 y se utilizará el año 2017 para validar el modelo:
Utilizamos la función lm():
modelo_lm<-lm(und~dia_semana+mes+Feriado,data=datos_ventas_autos,subset = (Fecha<="2016-12-31"))
Veamos el resultado:
summary(modelo_lm)
##
## Call:
## lm(formula = und ~ dia_semana + mes + Feriado, data = datos_ventas_autos,
## subset = (Fecha <= "2016-12-31"))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1077.39 -116.40 -20.32 108.56 2096.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 682.415 21.479 31.771 < 2e-16 ***
## dia_semana.L -698.583 16.790 -41.608 < 2e-16 ***
## dia_semana.Q -28.194 16.674 -1.691 0.091025 .
## dia_semana.C 729.146 16.603 43.915 < 2e-16 ***
## dia_semana^4 491.497 16.509 29.772 < 2e-16 ***
## dia_semana^5 8.992 16.495 0.545 0.585701
## dia_semana^6 -225.718 16.492 -13.687 < 2e-16 ***
## mesFebruary 117.821 31.003 3.800 0.000149 ***
## mesMarch 174.696 30.267 5.772 9.21e-09 ***
## mesApril 168.802 30.511 5.533 3.62e-08 ***
## mesMay 147.077 30.258 4.861 1.27e-06 ***
## mesJune 145.742 30.511 4.777 1.93e-06 ***
## mesJuly 153.337 30.261 5.067 4.45e-07 ***
## mesAugust 167.486 30.258 5.535 3.56e-08 ***
## mesSeptember 165.400 30.570 5.411 7.12e-08 ***
## mesOctober 150.525 30.273 4.972 7.25e-07 ***
## mesNovember 171.269 30.509 5.614 2.29e-08 ***
## mesDecember 464.365 30.259 15.346 < 2e-16 ***
## Feriado -819.039 28.770 -28.469 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 1808 degrees of freedom
## Multiple R-squared: 0.7696, Adjusted R-squared: 0.7673
## F-statistic: 335.5 on 18 and 1808 DF, p-value: < 2.2e-16
y0_tr<-mean(datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"])
r0_tr<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y0_tr
R0_tr<-mean(r0_tr^2)
y_pred_tr_lm<-predict(modelo_lm)
r_tr_lm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_lm
R_tr_lm<-mean(r_tr_lm^2)
R2_lm<-1-R_tr_lm/R0_tr
print(R2_lm)
## [1] 0.7695787
modelo_glm<-glm(und~dia_semana+mes+Feriado,data=datos_ventas_autos,subset = (Fecha<="2016-12-31"),family = "poisson")
Veamos el resultado:
summary(modelo_glm)
##
## Call:
## glm(formula = und ~ dia_semana + mes + Feriado, family = "poisson",
## data = datos_ventas_autos, subset = (Fecha <= "2016-12-31"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -42.581 -3.415 -1.120 2.933 45.878
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.795637 0.005850 990.71 <2e-16 ***
## dia_semana.L -2.559408 0.013146 -194.69 <2e-16 ***
## dia_semana.Q 0.101049 0.002206 45.81 <2e-16 ***
## dia_semana.C 2.800453 0.014192 197.32 <2e-16 ***
## dia_semana^4 3.120679 0.019455 160.40 <2e-16 ***
## dia_semana^5 1.858125 0.015167 122.51 <2e-16 ***
## dia_semana^6 0.524532 0.007138 73.49 <2e-16 ***
## mesFebruary 0.168204 0.004353 38.64 <2e-16 ***
## mesMarch 0.219482 0.004304 51.00 <2e-16 ***
## mesApril 0.217067 0.004310 50.36 <2e-16 ***
## mesMay 0.208480 0.004292 48.57 <2e-16 ***
## mesJune 0.203179 0.004353 46.67 <2e-16 ***
## mesJuly 0.207886 0.004270 48.68 <2e-16 ***
## mesAugust 0.234499 0.004276 54.84 <2e-16 ***
## mesSeptember 0.226172 0.004252 53.19 <2e-16 ***
## mesOctober 0.208291 0.004251 48.99 <2e-16 ***
## mesNovember 0.235782 0.004310 54.70 <2e-16 ***
## mesDecember 0.545568 0.004011 136.03 <2e-16 ***
## Feriado -5.662272 0.059671 -94.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 913454 on 1826 degrees of freedom
## Residual deviance: 108666 on 1808 degrees of freedom
## AIC: 121839
##
## Number of Fisher Scoring iterations: 8
Obtengamos el (pseudo) \(R^2\) para el modelo Poisson:
y_pred_tr_glm<-predict(modelo_glm,type="response")
r_tr_glm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_glm
R_tr_glm<-mean(r_tr_glm^2)
R2_tr_glm<-1-R_tr_glm/R0_tr
Veamos el resultado:
print(R2_tr_glm)
## [1] 0.807181
Creación de un dataframe para usar con plotly:
resultados_lm_glm<-data.frame(Fecha= datos_ventas_autos$Fecha[datos_ventas_autos$Fecha<="2016-12-31"], und=datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"],
pred_lm=y_pred_tr_lm,
pred_glm=y_pred_tr_glm,
res_lm=r_tr_lm,
res_glm=r_tr_glm)
plot_ly (data=resultados_lm_glm,
x = ~Fecha,
y = ~und,
type = "scatter" ,mode = "lines",
name='Real',
line=list(width=1,color='rgb(205, 12, 24)'))%>%
add_trace(y= ~pred_lm,
name='Modelo lineal general',
line=list(width=1,color='rgb(22, 96, 167)'))%>%
add_trace(y= ~pred_glm,
name='Modelo Poisson',
line=list(width=1,color='rgb(255, 51, 0)'))%>%
layout(title='Registros RUNT',
xaxis=list(title="Fecha"),
yaxis=list(title="Unidades"),
legend = list(x = 0.75, y = 0.9))
plot_ly (data=resultados_lm_glm,
x = ~und,
y = ~pred_lm,
text = ~Fecha,
type = "scatter" ,mode="markers",
name='Modelo lineal general',
marker=list(size=3,color='rgb(22, 96, 167)'))%>%
add_trace(y= ~pred_glm,
text = ~Fecha,
name='Modelo lineal Poisson',
marker=list(size=3,color='rgb(255, 51, 0)'))%>%
add_trace(x=c(-500:3500),y=c(-500:3500),
mode="lines",text=rep(NA,4001),
name="Identidad")%>%
layout(title='Registros RUNT',
xaxis=list(title="Observados"),
yaxis=list(title="Predichos"),
legend = list(x = 0.75, y = 0.9))
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
plot_ly (data=resultados_lm_glm,
x = ~und,
y = ~res_lm,
type = "scatter" ,mode="markers",
name='Modelo lineal general',
text = ~Fecha,
marker=list(size=3,color='rgb(22, 96, 167)'))%>%
add_trace(y= ~res_glm,
name='Modelo Poisson',
text = ~Fecha,
marker=list(size=3,color='rgb(255, 51, 0)'))%>%
layout(title='Registros RUNT',
xaxis=list(title="Observados"),
yaxis=list(title="Residuales"),
legend = list(x = 0.75, y = 0.9))
datos_val<-subset(datos_ventas_autos,subset=(Fecha>="2017-01-01" & Fecha<="2017-12-31"))
y_pred_vl_lm<-predict(modelo_lm,newdata = datos_val)
y_pred_vl_glm<-predict(modelo_glm,type="response",newdata = datos_val)
r_vl_lm<-datos_val$und-y_pred_vl_lm
r_vl_glm<-datos_val$und-y_pred_vl_glm
r0_vl<-datos_val$und-y0_tr
R_vl_lm<-mean(r_vl_lm^2)
R_vl_glm<-mean(r_vl_glm^2)
R0_vl<-mean(r0_vl^2)
R2_vl_lm<-1-R_vl_lm/R0_vl
R2_vl_glm<-1-R_vl_glm/R0_vl
print(R2_vl_lm)
## [1] 0.5184184
print(R2_vl_glm)
## [1] 0.5511457
resultados_lm_glm_val<-data.frame(Fecha=datos_val$Fecha, und=datos_val$und,
pred_lm=y_pred_vl_lm,
pred_glm=y_pred_vl_glm,
res_lm=r_vl_lm,
res_glm=r_vl_glm)
plot_ly (data=resultados_lm_glm_val,
x = ~Fecha,
y = ~und,
type = "scatter" ,mode = "lines",
name='Real',
line=list(width=1,color='rgb(205, 12, 24)'))%>%
add_trace(y= ~pred_lm,
name='Modelo lineal general',
line=list(width=1,color='rgb(22, 96, 167)'))%>%
add_trace(y= ~pred_glm,
name='Modelo Poisson',
line=list(width=1,color='rgb(255, 51, 0)'))%>%
layout(title='Registros RUNT (Validación)',
xaxis=list(title="Fecha"),
yaxis=list(title="Unidades"),
legend = list(x = 0.75, y = 0.9))
Creación de un vector que indica si la fecha es el último día hábil del mes. Los días hábiles son aquellos que no sean ni sábado ni domingo ni festivo:
dias_habiles<-ifelse(!(datos_ventas_autos$dia_semana %in% c("Saturday","Monday")) |
datos_ventas_autos$Feriado==0,1,0)
con_mod<-datos_ventas_autos$Consecutivo*dias_habiles
ultimo_dia_habil<-aggregate(con_mod~datos_ventas_autos$year*datos_ventas_autos$mes,FUN=max)
datos_ventas_autos$ultimo_dia_habil<-(datos_ventas_autos$Consecutivo %in% ultimo_dia_habil$con_mod)*1
Ahora incluyamos esta variable en el modelo:
modelo_glm<-glm(und~dia_semana+mes+Feriado+ultimo_dia_habil,
data=datos_ventas_autos,
subset = (Fecha<="2016-12-31"),
family = "poisson")
summary(modelo_glm)
##
## Call:
## glm(formula = und ~ dia_semana + mes + Feriado + ultimo_dia_habil,
## family = "poisson", data = datos_ventas_autos, subset = (Fecha <=
## "2016-12-31"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -51.846 -3.189 -0.814 3.231 45.198
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.782265 0.005853 987.99 <2e-16 ***
## dia_semana.L -2.561141 0.013146 -194.83 <2e-16 ***
## dia_semana.Q 0.100992 0.002206 45.78 <2e-16 ***
## dia_semana.C 2.799723 0.014192 197.27 <2e-16 ***
## dia_semana^4 3.119964 0.019455 160.36 <2e-16 ***
## dia_semana^5 1.857067 0.015167 122.44 <2e-16 ***
## dia_semana^6 0.521896 0.007138 73.12 <2e-16 ***
## mesFebruary 0.164640 0.004354 37.81 <2e-16 ***
## mesMarch 0.220647 0.004304 51.27 <2e-16 ***
## mesApril 0.214181 0.004311 49.69 <2e-16 ***
## mesMay 0.208480 0.004292 48.57 <2e-16 ***
## mesJune 0.206836 0.004354 47.51 <2e-16 ***
## mesJuly 0.205492 0.004270 48.12 <2e-16 ***
## mesAugust 0.235209 0.004276 55.01 <2e-16 ***
## mesSeptember 0.225107 0.004252 52.94 <2e-16 ***
## mesOctober 0.205699 0.004252 48.38 <2e-16 ***
## mesNovember 0.235935 0.004310 54.74 <2e-16 ***
## mesDecember 0.543046 0.004011 135.39 <2e-16 ***
## Feriado -5.647972 0.059671 -94.65 <2e-16 ***
## ultimo_dia_habil 0.349682 0.003884 90.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 913454 on 1826 degrees of freedom
## Residual deviance: 101355 on 1807 degrees of freedom
## AIC: 114531
##
## Number of Fisher Scoring iterations: 8
Veamos cómo cambia el \(R^2\):
y_pred_tr_glm<-predict(modelo_glm,type="response")
r_tr_glm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_glm
R_tr_glm<-mean(r_tr_glm^2)
R2_tr_glm<-1-R_tr_glm/R0_tr
print(R2_tr_glm)
## [1] 0.8183819
Y en validación:
datos_val<-subset(datos_ventas_autos,subset=(Fecha>="2017-01-01" & Fecha<="2017-12-31"))
y_pred_vl_glm<-predict(modelo_glm,type="response",newdata = datos_val)
r_vl_glm<-datos_val$und-y_pred_vl_glm
R_vl_glm<-mean(r_vl_glm^2)
R2_vl_glm<-1-R_vl_glm/R0_vl
print(R2_vl_glm)
## [1] 0.6444367